#------ Check the Appedix of Kim et al. (2015) for the model specification and hyper-parameters.

data {
    int<lower=0> n;
    int<lower=0> cov;
    real out[n];
    row_vector[cov] x[n];
    int C;
    real<lower=0> precision;
    real mu_sub_mean;
}
parameters {
    real gamma0[C];
    vector[cov] gamma1;
    vector<lower=0, upper=1>[C] V;
    real<lower=0> W[C];
    real<lower=0.3> alpha;
    real<lower=0> w_a;
    real mu_sub;
    real<lower=0> w_sub;
    real<lower=0> w_suba;
}
transformed parameters {
    real<lower=0> w_b;
    real<lower=0> w_subb;
    vector<lower=0, upper=1>[C] psi;
    w_b <- w_a / (2*precision);
    w_subb <- w_suba / (2*precision);
    psi[1] <- V[1];
    for (j in 2:C){
        psi[j] <- V[j] * (1 - V[j - 1]) * psi[j - 1] / V[j - 1];
    }
    for (j in 1:C){
        psi[j] <- psi[j] / sum(psi);
    }
}
model {
    real ps[C];
    alpha ~ gamma(0.1, 0.1);
    for (j in 1:C){
        V[j] ~ beta(1, alpha);
    }
    w_a ~ uniform(1,10);
    w_suba ~ uniform(1,10);
    w_sub ~ gamma(w_suba, w_subb);
    mu_sub ~ normal(mu_sub_mean, 1/(2*precision));
    for (j in 1:cov){
        gamma1[j] ~ normal(0, 100);
    }
    for (j in 1:C){
        gamma0[j] ~ normal(mu_sub, 1/w_sub);
        W[j] ~ gamma(w_a, w_b);
    }
    for (j in 1:n){
        for (k in 1:C){
            ps[k] <- log(psi[k]) + normal_log(out[j], gamma0[k]+x[j]*gamma1, 1/W[k]);
        }
        increment_log_prob(log_sum_exp(ps));
    }
  
}
